rm(list = ls())
###

###
library(tidyverse) # 1.3.1
library(stargazer) # 5.2.2
library(scales) # 1.1.1
library(modelsummary) # 0.7.0
library(vtable) # 1.3.1
library(gtsummary) # 1.4.1
library(gridExtra) # 2.3
library(jtools) # 2.1.3

# Load data
data <- read.csv("perspective_taking_merged.csv")

## Filter out those who refused informed consent 
# + failed attention checks

data %>% 
  filter(informed_consent == 1 & attention_check1 == 3) %>%
  mutate_all(na_if,"") %>%
  filter(age >= 18) %>%
  filter(!is.na(attention_check1)) -> data



#### Use favorability pre and post
data %>%
  mutate(prejudice_pre = recode(undoc_immigrants, # higher value = more favorable
                              `0` = 0,
                              `1` = .25,
                              `2` = .50,
                              `3` = .75,
                              `4` = 1),
         prejudice_post = recode(post_s1_prejudice, ## flip, lower values => less upset = more favorable
                                 `0` = 1, 
                                 `1` = .8, 
                                 `2` = .6,
                                 `3` = .4,
                                 `4` = .2,
                                 `5` = 0)) %>%
  mutate(treatment = factor(treatment,
                              levels = c("Control Undoc",
                                         "Treat Undoc",
                                         "Cancer",
                                         "Immigrant",
                                         "Control Low Inc",
                                         "Treat Low Inc")) ) -> pre_post

## Calculate mean favorability pre and post
pre_post %>%
  group_by(treatment) %>%
  summarise(n = n(),
            mean_pre = mean(prejudice_pre, na.rm = T),
            se_pre = sd(prejudice_pre, na.rm = T) / sqrt(n()),
            mean_post = mean(prejudice_post, na.rm = T),
            se_post = sd(prejudice_post, na.rm=T) / sqrt(n())) %>%
  filter(!is.na(treatment))

##### Visualize data
## Get means for pre-treatment
pre_post %>%
  group_by(treatment) %>%
  summarise(n = n(),
            mean = mean(prejudice_pre, na.rm = T),
            se = sd(prejudice_pre, na.rm = T) / sqrt(n())) %>%
  filter(!is.na(treatment) ) %>%
  mutate(indicator = "pre") -> pre

# Get means for post-treatment
pre_post %>%
  group_by(treatment) %>%
  summarise(n = n(),
            mean = mean(prejudice_post, na.rm = T),
            se = sd(prejudice_post, na.rm=T) / sqrt(n())) %>%
  filter(!is.na(treatment) ) %>%
  mutate(indicator = "post") -> post

# Bind data
pre_post_means <- as.data.frame(rbind(pre, post))

####################################      
###### Figure A16 ################## 
####################################      
pre_post_means %>%
  mutate(order = recode(treatment, ## order the treatments
                        "Control Undoc" = 1,
                        "Treat Undoc" = 2,
                        "Cancer" = 3,
                        "Immigrant" = 4,
                        "Control Low Inc" = 5,
                        "Treat Low Inc" = 6),
         indicator = factor(indicator, # create indicator
                            levels = c("pre",
                                       "post")),
         treated = factor(case_when(treatment == "Control Undoc" |        # label treatment/control
                               treatment == "Control Low Inc" ~ "Control",
                             TRUE ~ "Treated"),
                          levels = c("Control",
                                     "Treated")) ) %>%
  ggplot() +
  geom_pointrange(aes(reorder(treatment, -order), 
                      y = mean,
                      ymin = mean - 2*se,
                      ymax = mean + 2*se,
                      color = as.factor(treated),
                      shape = as.factor(treated)),
                  position = position_dodge2(w = .5)) +
  coord_flip() +
  theme_bw() +
  facet_wrap(. ~ indicator,
             labeller = as_labeller(c(`pre` = "Pre-Treatment",
                                      `post` = "Post-Treatment"))
             ) +
  theme(legend.position = "top") +
  scale_y_continuous(labels = scales::percent) +
  scale_x_discrete(labels = c("Control Undoc" = "Control\n+Undocumented",
                              "Treat Undoc" = "Treatment\n+Undocumented\n+COVID-19",
                              "Cancer" = "Treatment\n+Undocumented\n+Cancer",
                              "Immigrant" = "Treatment\n+Immigrant\n+COVID-19",
                              "Control Low Inc" = "Control\n+Low Income",
                              "Treat Low Inc" = "Treatment\n+Low Income\n+COVID-19") ) +
  scale_shape_discrete(name = "",
                       labels = c("Control",
                                  "Treatment")) +
  scale_color_manual(name = "",
                     labels = c("Control",
                                "Treatment"),
                     values = c("dark blue",
                                "red3")) +
  labs(x = "", y = "Mean Favorability Toward Undoc Imm") 
ggsave("appendix_figure_a16_final.jpeg")

####################################      
###### Figure A17 ################## 
####################################

# create figure
plot_summs(list(
  "Unadjusted" = lm(prejudice_post ~ treatment, pre_post),
  "Adjusted" = lm(prejudice_post ~ treatment
                  + age + gender + income + educ + ethnicity + pid_3level, pre_post)),
  coefs = c("Treatment\n+Undocumented\n+COVID-19" = "treatmentTreat Undoc",
            "Treatment\n+Undocumented\n+Cancer" = "treatmentCancer",
            "Treatment\n+Immigrant\n+COVID-19" = "treatmentImmigrant",
            "Control\n+Low Income" = "treatmentControl Low Inc",
            "Treatment\n+Low Income\n+COVID-19" = "treatmentTreat Low Inc"),
  inner_ci_level = .90, robust = T,
  colors = c("dark blue", "light blue")
) +
  theme(legend.position = "top")
ggsave("appendix_figure_a17_final.jpeg")


